#!/usr/bin/env julia

"""
    test_observation_processing.jl

Comprehensive test suite for observation processing components in GSI.
Tests all major observation types, processing pipelines, and integration
with core GSI analysis algorithms.

Test Coverage:
- Observation type construction and validation
- Basic processing operations (interpolation, quality control)
- Integration with forward operators and bias correction
- Memory management and performance characteristics
- Error handling and edge cases
"""

using Test
using GSICoreAnalysis
using LinearAlgebra
using Random
using Statistics
using Dates

# Import observation processing modules
using GSICoreAnalysis.ObservationProcessing

"""Test observation type creation and basic properties"""
@testset "Observation Types Construction" begin
    println("  Testing observation type construction...")
    
    # Test conventional observation creation
    @testset "Conventional Observations" begin
        # Surface pressure observation
        ps_obs = ConventionalObservation{Float64}(
            obs_type = "surface_pressure",
            value = 1013.25,
            error = 2.0,
            latitude = 40.0,
            longitude = -74.0,
            elevation = 10.0,
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        @test ps_obs.obs_type == "surface_pressure"
        @test ps_obs.value == 1013.25
        @test ps_obs.error == 2.0
        @test ps_obs.latitude == 40.0
        @test ps_obs.longitude == -74.0
        @test ps_obs.elevation == 10.0
        @test ps_obs.qc_flag == 0  # Default should be good
        
        # Temperature observation
        temp_obs = ConventionalObservation{Float64}(
            obs_type = "temperature",
            value = 285.5,
            error = 1.5,
            latitude = 45.5,
            longitude = -122.7,
            elevation = 500.0,
            pressure = 850.0,
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        @test temp_obs.obs_type == "temperature"
        @test temp_obs.pressure == 850.0
        @test !isnothing(temp_obs.pressure)
    end
    
    # Test satellite observation creation
    @testset "Satellite Observations" begin
        # AMSU-A radiance observation
        amsu_obs = SatelliteObservation{Float64}(
            obs_type = "amsu_a",
            channel = 5,
            value = 245.2,
            error = 1.2,
            latitude = 35.0,
            longitude = 139.0,
            satellite_zenith_angle = 15.3,
            satellite_azimuth_angle = 180.2,
            scan_position = 15,
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        @test amsu_obs.obs_type == "amsu_a"
        @test amsu_obs.channel == 5
        @test amsu_obs.satellite_zenith_angle == 15.3
        @test amsu_obs.scan_position == 15
        
        # GPS RO observation
        gps_obs = SatelliteObservation{Float64}(
            obs_type = "gps_ro",
            channel = 1,
            value = 1.25e-4,  # Refractivity
            error = 5.0e-6,
            latitude = 60.0,
            longitude = 25.0,
            impact_parameter = 6371.0e3 + 15.0e3,  # 15 km altitude
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        @test gps_obs.obs_type == "gps_ro"
        @test gps_obs.impact_parameter ≈ 6.386e6
    end
    
    # Test observation collection
    @testset "Observation Collections" begin
        obs_list = ObservationData{Float64}()
        
        # Add conventional observations
        ps_obs = ConventionalObservation{Float64}(
            obs_type = "surface_pressure",
            value = 1015.0,
            error = 2.0,
            latitude = 30.0,
            longitude = -90.0,
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        temp_obs = ConventionalObservation{Float64}(
            obs_type = "temperature", 
            value = 288.0,
            error = 1.5,
            latitude = 30.0,
            longitude = -90.0,
            pressure = 900.0,
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        add_observation!(obs_list, ps_obs)
        add_observation!(obs_list, temp_obs)
        
        @test length(obs_list.conventional) == 2
        @test length(obs_list.satellite) == 0
        
        # Test observation counts by type
        counts = count_observations_by_type(obs_list)
        @test counts["surface_pressure"] == 1
        @test counts["temperature"] == 1
    end
end

"""Test observation processing operations"""
@testset "Basic Processing Operations" begin
    println("  Testing basic processing operations...")
    
    @testset "Spatial Operations" begin
        # Create test observation
        obs = ConventionalObservation{Float64}(
            obs_type = "temperature",
            value = 290.0,
            error = 1.0,
            latitude = 40.0,
            longitude = -105.0,
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        # Test distance calculations
        lat2, lon2 = 41.0, -104.0
        distance = great_circle_distance(obs.latitude, obs.longitude, lat2, lon2)
        @test distance > 0.0
        @test distance < 200.0e3  # Should be less than 200 km
        
        # Test within radius check
        @test is_within_radius(obs, 40.5, -104.5, 100.0e3)  # Within 100 km
        @test !is_within_radius(obs, 45.0, -110.0, 200.0e3)  # More than 200 km away
    end
    
    @testset "Temporal Operations" begin
        base_time = DateTime(2024, 1, 1, 12, 0, 0)
        
        obs1 = ConventionalObservation{Float64}(
            obs_type = "temperature",
            value = 285.0,
            error = 1.0,
            latitude = 40.0,
            longitude = -100.0,
            time = base_time
        )
        
        obs2 = ConventionalObservation{Float64}(
            obs_type = "temperature", 
            value = 287.0,
            error = 1.0,
            latitude = 40.0,
            longitude = -100.0,
            time = base_time + Minute(30)
        )
        
        # Test time differences
        time_diff = temporal_distance(obs1, obs2)
        @test time_diff == 1800.0  # 30 minutes in seconds
        
        # Test within time window
        @test is_within_time_window(obs1, base_time, 3600.0)  # Within 1 hour
        @test is_within_time_window(obs2, base_time, 3600.0)  # Within 1 hour
        @test !is_within_time_window(obs2, base_time, 1500.0) # Outside 25 minutes
    end
end

"""Test quality control algorithms"""
@testset "Quality Control" begin
    println("  Testing quality control algorithms...")
    
    @testset "Background Check" begin
        # Create test observation
        obs = ConventionalObservation{Float64}(
            obs_type = "temperature",
            value = 290.0,
            error = 1.5,
            latitude = 40.0,
            longitude = -100.0,
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        # Test background departure check
        background_value = 288.5
        departure = obs.value - background_value
        
        # Should pass QC with small departure
        @test abs(departure) < 5.0 * obs.error
        
        # Test gross error detection
        obs_bad = ConventionalObservation{Float64}(
            obs_type = "temperature",
            value = 350.0,  # Unrealistic temperature
            error = 1.5,
            latitude = 40.0,
            longitude = -100.0,
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        bad_departure = obs_bad.value - background_value
        @test abs(bad_departure) > 10.0 * obs_bad.error  # Should be flagged
    end
    
    @testset "Buddy Check" begin
        # Create nearby observations for buddy check
        obs_list = ObservationData{Float64}()
        
        # Good central observation
        central_obs = ConventionalObservation{Float64}(
            obs_type = "temperature",
            value = 288.0,
            error = 1.0,
            latitude = 40.0,
            longitude = -100.0,
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        # Nearby supporting observations
        for i in 1:4
            lat_offset = (i-2.5) * 0.1  # ±0.15 degrees
            lon_offset = (i-2.5) * 0.1
            buddy_value = 288.0 + randn() * 0.5  # Small random variations
            
            buddy_obs = ConventionalObservation{Float64}(
                obs_type = "temperature",
                value = buddy_value,
                error = 1.0,
                latitude = 40.0 + lat_offset,
                longitude = -100.0 + lon_offset,
                time = DateTime(2024, 1, 1, 12, 0, 0)
            )
            
            add_observation!(obs_list, buddy_obs)
        end
        
        add_observation!(obs_list, central_obs)
        
        # Find buddies within radius
        buddies = find_buddy_observations(central_obs, obs_list.conventional, 50.0e3)
        @test length(buddies) >= 3  # Should find most of the nearby obs
        
        # Compute buddy statistics
        buddy_values = [obs.value for obs in buddies]
        buddy_mean = mean(buddy_values)
        buddy_std = std(buddy_values)
        
        @test abs(central_obs.value - buddy_mean) < 2.0 * buddy_std
    end
end

"""Test forward operator integration"""
@testset "Forward Operator Integration" begin
    println("  Testing forward operator integration...")
    
    @testset "Interpolation Operations" begin
        # Create a simple test grid
        config = AnalysisConfig{Float64}(
            grid_size = (8, 6, 4),
            ensemble_size = 0
        )
        
        sv = StateVector(config)
        
        # Initialize with a simple temperature field (linear gradient)
        fill!(sv.values, 0.0)
        temp_slice = get_variable_slice(sv, "temperature")
        
        for k in 1:config.grid_size[3]
            for j in 1:config.grid_size[2]
                for i in 1:config.grid_size[1]
                    # Simple linear temperature field: T = 300 - 20*(k-1)/3
                    idx = linear_index(i, j, k, config.grid_size)
                    temp_slice[idx] = 300.0 - 20.0 * (k-1) / (config.grid_size[3]-1)
                end
            end
        end
        
        # Test observation
        obs = ConventionalObservation{Float64}(
            obs_type = "temperature",
            value = 285.0,
            error = 1.5,
            latitude = 40.0,  # Middle of domain
            longitude = -100.0,  # Middle of domain  
            pressure = 850.0,  # Mid-level
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        # Test forward operator (H[x])
        h_x = apply_forward_operator(obs, sv, config)
        @test h_x > 280.0
        @test h_x < 300.0  # Should be within reasonable range
        
        # Test that interpolation is working
        @test !isnan(h_x)
        @test isfinite(h_x)
    end
    
    @testset "Tangent Linear and Adjoint" begin
        config = AnalysisConfig{Float64}(
            grid_size = (6, 6, 3),
            ensemble_size = 0
        )
        
        sv = StateVector(config)
        cv = ControlVector(config)
        
        # Initialize state and control vectors
        randn!(sv.values)
        randn!(cv.values)
        
        # Create test observation
        obs = ConventionalObservation{Float64}(
            obs_type = "temperature",
            value = 290.0,
            error = 2.0,
            latitude = 45.0,
            longitude = -110.0,
            pressure = 700.0,
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        # Test tangent linear operator H⋅dx
        h_tl = apply_tangent_linear_operator(obs, cv, sv, config)
        @test isfinite(h_tl)
        
        # Test adjoint operator H^T⋅y
        innovation = 5.0  # obs - H(x)
        apply_adjoint_operator!(obs, innovation, cv, sv, config)
        
        # Check that adjoint modified control vector
        cv_norm_after = norm_cv(cv)
        @test cv_norm_after > 0.0
        
        # Test adjoint property: <H⋅dx, y> = <dx, H^T⋅y>
        # This would require more sophisticated testing framework
    end
end

"""Test bias correction algorithms"""
@testset "Bias Correction" begin
    println("  Testing bias correction algorithms...")
    
    @testset "Satellite Radiance Bias Correction" begin
        # Create satellite observation
        sat_obs = SatelliteObservation{Float64}(
            obs_type = "amsu_a",
            channel = 6,
            value = 242.5,
            error = 1.8,
            latitude = 45.0,
            longitude = -120.0,
            satellite_zenith_angle = 25.0,
            scan_position = 20,
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        # Initialize bias correction parameters
        bias_params = SatelliteBiasCorrection{Float64}(
            satellite_id = "noaa-18",
            sensor_id = "amsu-a",
            channel = 6,
            scan_angle_coeffs = [0.1, -0.05, 0.02],
            airmass_coeffs = [0.2, -0.1],
            constant_bias = -0.5
        )
        
        # Apply bias correction
        corrected_value = apply_bias_correction(sat_obs, bias_params)
        @test corrected_value != sat_obs.value  # Should be different
        @test abs(corrected_value - sat_obs.value) < 5.0  # But not drastically different
        
        # Test that bias correction is reasonable
        @test corrected_value > 200.0  # Still a reasonable radiance value
        @test corrected_value < 300.0
    end
    
    @testset "Conventional Observation Bias Correction" begin  
        # Create surface pressure observation
        conv_obs = ConventionalObservation{Float64}(
            obs_type = "surface_pressure",
            value = 1015.3,
            error = 1.5,
            latitude = 40.0,
            longitude = -95.0,
            elevation = 300.0,
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        # Simple elevation-based bias correction
        elevation_bias = -0.12 * conv_obs.elevation / 1000.0  # -0.12 hPa per km
        corrected_value = conv_obs.value + elevation_bias
        
        @test corrected_value < conv_obs.value  # Should reduce for positive elevation
        @test abs(corrected_value - conv_obs.value) < 1.0  # Small correction
    end
end

"""Test spatial processing (thinning and super-observations)"""
@testset "Spatial Processing" begin
    println("  Testing spatial processing...")
    
    @testset "Observation Thinning" begin
        # Create dense observation set
        obs_list = ObservationData{Float64}()
        
        # Add many observations in a small area
        base_lat, base_lon = 35.0, -78.0
        for i in 1:20
            for j in 1:20
                lat = base_lat + (i-10) * 0.01  # 0.2° x 0.2° box
                lon = base_lon + (j-10) * 0.01
                
                obs = ConventionalObservation{Float64}(
                    obs_type = "temperature",
                    value = 288.0 + randn() * 1.0,
                    error = 1.5,
                    latitude = lat,
                    longitude = lon,
                    time = DateTime(2024, 1, 1, 12, 0, 0)
                )
                
                add_observation!(obs_list, obs)
            end
        end
        
        @test length(obs_list.conventional) == 400
        
        # Apply thinning with 5 km resolution
        thinned_obs = thin_observations(obs_list.conventional, 5.0e3)
        
        # Should significantly reduce observation count
        @test length(thinned_obs) < 50
        @test length(thinned_obs) > 0
        
        # Verify minimum distances
        for i in 1:length(thinned_obs)-1
            for j in i+1:length(thinned_obs)
                dist = great_circle_distance(
                    thinned_obs[i].latitude, thinned_obs[i].longitude,
                    thinned_obs[j].latitude, thinned_obs[j].longitude
                )
                @test dist >= 4.5e3  # Allow small tolerance
            end
        end
    end
    
    @testset "Super-Observation Creation" begin
        # Create clustered observations for super-obs
        obs_cluster = ConventionalObservation{Float64}[]
        
        center_lat, center_lon = 40.0, -100.0
        true_temp = 285.0
        
        # Add observations with small spatial spread
        for i in 1:10
            lat_noise = randn() * 0.05  # ~5 km standard deviation
            lon_noise = randn() * 0.05
            temp_noise = randn() * 0.5   # 0.5 K measurement noise
            
            obs = ConventionalObservation{Float64}(
                obs_type = "temperature",
                value = true_temp + temp_noise,
                error = 1.0,
                latitude = center_lat + lat_noise,
                longitude = center_lon + lon_noise, 
                time = DateTime(2024, 1, 1, 12, 0, 0)
            )
            
            push!(obs_cluster, obs)
        end
        
        # Create super-observation
        super_obs = create_superobservation(obs_cluster)
        
        @test super_obs.obs_type == "temperature"
        @test abs(super_obs.latitude - center_lat) < 0.1
        @test abs(super_obs.longitude - center_lon) < 0.1
        @test abs(super_obs.value - true_temp) < 1.0  # Should be close to truth
        @test super_obs.error < 1.0  # Should have reduced error
        @test super_obs.error > 0.1  # But not unrealistically small
    end
end

"""Test memory management and performance"""
@testset "Memory Management" begin
    println("  Testing memory management...")
    
    @testset "Large Observation Sets" begin
        # Create large observation dataset
        n_obs = 10_000
        obs_list = ObservationData{Float64}()
        
        # Add many observations
        for i in 1:n_obs
            lat = -90.0 + 180.0 * rand()
            lon = -180.0 + 360.0 * rand()
            
            obs = ConventionalObservation{Float64}(
                obs_type = "temperature",
                value = 250.0 + 50.0 * rand(),
                error = 0.5 + 1.5 * rand(),
                latitude = lat,
                longitude = lon,
                time = DateTime(2024, 1, 1, 12, 0, 0)
            )
            
            add_observation!(obs_list, obs)
        end
        
        @test length(obs_list.conventional) == n_obs
        
        # Test memory cleanup
        clear_observations!(obs_list)
        @test length(obs_list.conventional) == 0
        @test length(obs_list.satellite) == 0
    end
    
    @testset "Processing Performance" begin
        # Create moderate-size dataset for performance testing
        n_obs = 1000
        obs_list = ConventionalObservation{Float64}[]
        
        for i in 1:n_obs
            obs = ConventionalObservation{Float64}(
                obs_type = "temperature",
                value = 288.0 + randn() * 5.0,
                error = 1.0 + 0.5 * rand(),
                latitude = -45.0 + 90.0 * rand(),
                longitude = -180.0 + 360.0 * rand(),
                time = DateTime(2024, 1, 1, 12, 0, 0)
            )
            
            push!(obs_list, obs)
        end
        
        # Time the thinning operation
        start_time = time()
        thinned_obs = thin_observations(obs_list, 50.0e3)  # 50 km resolution
        processing_time = time() - start_time
        
        @test processing_time < 5.0  # Should complete within 5 seconds
        @test length(thinned_obs) > 0
        @test length(thinned_obs) <= length(obs_list)
    end
end

"""Test error handling and edge cases"""  
@testset "Error Handling" begin
    println("  Testing error handling and edge cases...")
    
    @testset "Invalid Observation Data" begin
        # Test invalid coordinates
        @test_throws AssertionError ConventionalObservation{Float64}(
            obs_type = "temperature",
            value = 288.0,
            error = 1.0,
            latitude = 95.0,  # Invalid latitude
            longitude = -100.0,
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        @test_throws AssertionError ConventionalObservation{Float64}(
            obs_type = "temperature",
            value = 288.0,
            error = 1.0,
            latitude = 45.0,
            longitude = 185.0,  # Invalid longitude
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        # Test negative error
        @test_throws AssertionError ConventionalObservation{Float64}(
            obs_type = "temperature", 
            value = 288.0,
            error = -1.0,  # Invalid error
            latitude = 45.0,
            longitude = -100.0,
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
    end
    
    @testset "Empty Observation Lists" begin
        empty_list = ObservationData{Float64}()
        
        # Test operations on empty lists
        counts = count_observations_by_type(empty_list)
        @test length(counts) == 0
        
        # Test thinning empty list
        thinned = thin_observations(empty_list.conventional, 10.0e3)
        @test length(thinned) == 0
        
        # Test finding buddies with no observations
        test_obs = ConventionalObservation{Float64}(
            obs_type = "temperature",
            value = 288.0,
            error = 1.0,
            latitude = 40.0,
            longitude = -100.0,
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        buddies = find_buddy_observations(test_obs, empty_list.conventional, 50.0e3)
        @test length(buddies) == 0
    end
    
    @testset "Extreme Values" begin
        # Test observation with extreme but valid values
        extreme_obs = ConventionalObservation{Float64}(
            obs_type = "temperature",
            value = 183.0,  # Very cold but physically possible
            error = 0.1,    # Very precise
            latitude = -89.9,  # Near south pole
            longitude = 179.9,  # Near date line
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        @test extreme_obs.value == 183.0
        @test extreme_obs.latitude == -89.9
        
        # Test operations with extreme values
        normal_obs = ConventionalObservation{Float64}(
            obs_type = "temperature",
            value = 288.0,
            error = 1.0,
            latitude = 0.0,
            longitude = 0.0, 
            time = DateTime(2024, 1, 1, 12, 0, 0)
        )
        
        # Distance should be computable
        distance = great_circle_distance(
            extreme_obs.latitude, extreme_obs.longitude,
            normal_obs.latitude, normal_obs.longitude
        )
        @test distance > 0.0
        @test isfinite(distance)
    end
end

println("  ✓ All observation processing tests completed successfully")